Inverse approach to atomistic modeling: Applications to a-Si:H and g-GeSe2 
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We discuss an inverse approach for atomistic modeling of glassy materials. The focus is on 
structural modeling and electronic properties of hydrogenated amorphous silicon and glassy GeSe2 
alloy. The work is based upon a new approach "experimentally constrained molecular relaxation 
(ECMR)". Unlike conventional approaches (such as molecular dynamics (MD) and Monte Carlo 
simulations(MC), where a potential function is specified and the system evolves either determinis- 
tically (MD) or stochastically (MC), we develop a novel scheme to model structural configurations 
using experimental data in association with density functional calculations. We have applied this 
approach to model hydrogenated amorphous silicon and glassy GeSe2. The electronic and struc- 
tural properties of these models are compared with experimental data and models obtained from 
conventional molecular dynamics simulation. 

PACS numbers: 71.23.Cq, 71.15.Mb, 71.23.An 



INTRODUCTION 

In conventional electronic structure problem one starts 
with a formulation of a model that consists of a set of 
atomic or molecular coordinates and an interacting po- 
tential or Hamiltonian. The electronic density of states 
is obtained either by solving the Schrodinger equation 
within the first-principles density functional formalism 
or by constructing a semi-empirical Hamiltonian to com- 
pute the electronic eigenstates, total energy, and response 
functions to compare with experiments For a crys- 
talline system, such an approach is fairly straightforward 
(at least in principle), although in many cases it can 
be computationally very expensive. Disordered materi- 
als, however, pose the problem that is highly non-trivial 
due to lack of structural information about the materi- 
als. Amorphous materials are archetypal examples of this 
class, which are characterized by presence of short and 
medium range order ranging from few angstroms to sev- 
eral nanometers . Since the position of an atom is not 
known in the amorphous state, one needs to proceed by 
simultaneously optimizing the geometry and calculating 
the electronic properties and structure for that geome- 
try. While there exists a number of methods that can 
address this problem with a varying degree of accuracy 
(such as Car-Parrinello Q and first-principles molecu- 
lar dynamics almost all cases where the problem is 
characterized by long time and large length scale, these 
techniques are computationally overkill and largely in- 
applicable for problems that require realistic modeling 
of complex amorphous systems. The so-called order-N 
methods (that scale linearly with system size N) can deal 
with large number of atoms; however, their use is of- 
ten severely limited in structural modeling of amorphous 
materials due to lack of knowledge in spectral proper- 
ties (e.g. the nature of density matrix and the presence 



of spectral gap). First-principles molecular dynamics is, 
therefore, not practically viable for large scale modeling 
that are necessary to study some of the novel properties 
of materials associated with amorphous state (such as 
medium range order and diffusion in glasses). The use 
of empirical potentials make it possible to run molecular 
dynamics simulation for several thousands atoms with a 
time scale of the order of several nanoseconds, but for 
many complex systems empirical potentials are not reli- 
able enough to describe the geometry and local chemistry 
correctly. 

In this paper, we discuss applications of an alterna- 
tive approach to model materials using a combination of 
experimental data and a suitable force-field (either first- 
principles, tight-binding or empirical). Instead of taking 
a direct approach, we take an inverse approach where 
experimental data are enforced to build atomic config- 
urations that have the desired structural and electronic 
properties. Reverse Monte Carlo (RMC) is a classic ex- 
ample of such a method that has been discussed by sev- 
eral authors [4, iSi, i6|, iTj . Here we briefly mention an ex- 
tension of reverse Monte Carlo by merging experimental 
data using a suitable force-field known as Experimental 
Constrained Molecular Relaxation (ECMR), and discuss 
its application to hydrogenated amorphous silicon and 
glassy- GeSe2. 



EXPERIMENTALLY CONSTRAINED 
MOLECULAR RELAXATION (ECMR) 

The reverse Monte Carlo (RMC) method developed by 
McGreevy and colleagues is a classic example of an in- 
verse problem in materials modeling |^]. The method 
constructs a model configuration by making use of all 
available experimental information. The central idea is 
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to set up a generalised cost function containing as much 
information as possible, and then to minimize the func- 
tion for generating configurations toward exact agree- 
ment with one or more experiments. In order to reduce 
the solution space and to explore a limited part of phys- 
ical subspace of configuration, a number of chemical, ge- 
ometrical and topological constraints can be added. The 
mathematical structure of this problem is equivalent to 
the constraint optimization "traveling salesman" prob- 
lem: 
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The coefficients 77^- and Xj are the appropriate weight fac- 
tors for each data set (F) and the constraints (P). The 
quantity Q is the appropriate generalized variable asso- 
ciated with experimental data F{Q). In order to prevent 
the atoms getting too close to each other, a certain cut- 
off distance is imposed, which is typically of the order 
of interatomic spacing. RMC has been applied success- 
fully to variety of materials [4]. However, the principal 
difficulty is the lack of uniqueness of the method. In 
absence of information from high order correlation func- 
tions, the method produces a range of configurations that 
are consistent with input experimental data, but may not 
be physically meaningful. The imposition of topological 
or chemical constraints can ameliorate the problem, but 
cannot eliminate completely. 

The experimentally constrained molecular relaxation 
(ECMR) has been designed to overcome some of the 
problems above Q. Instead of relying on experimen- 
tal information and a set of constraints only, one addi- 
tionally employs an approximate energy functional to de- 
scribe the dynamics correctly by merging first-principles 
density functional (or semi-empirical) as well as experi- 
mental data. The purpose of the energy functional is to 
guide the system approximately in the augmented config- 
urational space defined by experimental data and other 
constraints. This largely reduces the number of unphys- 
ical solutions that are mathematically correct but fail to 
satisfy the dynamical behavior correctly. Introduction of 
an approximate energy functional constrains the system 
to evolve on a restricted hut more realistic energy sur- 
face, and thereby accelerates to converge toward reliable 
structural configurations during the course of simulation. 
The configuration obtained from the method is not only a 
minimum (metastable or global) of an appropriate energy 
functional but also consistent with the input experimen- 
tal information. Symbolically, we can write the grand 
penalty function in the ECMR as: 



S(Q,r)-e(Q) © 7£;(r) 
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that of the energy functional E{y) [9|. In the limit 7 
is infinitesimally small, the method reduces to a inverse 
method (RMC in the present case), whereas for a very 
large value of 7 the method is equivalent to a direct 
method of minimizing the total energy. 



APPLICATIONS: GLASSY-GESE2 AND A-SI:H 

Amorphous GeSe2 is a classic glass former and has in- 
teresting physical properties that are difficult to model 
via conventional molecular dynamics simulation [lo| . 
The material strongly shows the presence of intermediate 
range order in the form of a first sharp diffraction peak 
(FSDP) in Neutron diffraction measurements. The origin 
of this intermediate range order is generally attributed 
to the presence of tetrahedral motifs having edge- and 
corner-sharing topology. Raman spectroscopy and Neu- 
tron diffraction provide useful information about the 
topological structure of the material 0,[l3. The GeSea 
model simulated in our ECMR work consists of 647 
atoms. The static structure factor (Neutron-weighted) 
of the final model is plotted in figure [T] along with the 
experimental data and the data obtained by Cobb et 
al. in Ref. [l^. It is clear from the figure that experi- 
mental data fit with the theoretical values obtained from 
ECMR model reasonably good. The structure at large Q 
matches quite well, but at small Q there exists deviations 
from the experimental peak position. 

Since the network topology of g-GeSe2 is very much 
influenced by the presence of edge- and corner-sharing 
tetrahedra, we have studied the presence of such tetra- 
hedra in our model. As mentioned before, Raman and 
Neutron experiments 11, l3| revealed that about 33% to 



In equation ([2]), the symbol ® stands for direct sum 
of the configuration space of penalty function ^(Q) and 



40% of Ge atoms are involved in forming edge-sharing 
tetrahedra. The corresponding percentage in our model 
is found to be 38%. This is particularly important in 
view of the fact that such information is not included in 
our starting ECMR model. It is remarkable to note that 
imposition of partial pair correlation functions and first- 
principles relaxation via ECMR does introduce the cor- 
rect topological structure in the network. We have also 
observed that 81% of Ge atoms in our model are 4- fold 
coordinated of which approximately 75% form predomi- 
nant Ge-center motifs Ge{Se^)A while 6% are ethane-like 
Ge2 (5*61)6 units. The remaining Ge atoms are 3-fold 
coordinated and are mostly found to be bonded as Ge- 
Ses units. The experimental radial distributions (partial) 
provide the partial nearest neighbor coordination num- 
ber. The values for Ge-Ge, Ge-Se and Se-Se are given by 
0.25, 0.20 and 3.7 respectively, whereas the correspond- 
ing values from our ECMR model are 0.17, 0.30 and 3.68 
respectively. In addition to studying the structural prop- 
erties, we have also computed the electronic properties 
of the model. The electronic density of states provide an 
additional check and can be compared with X-ray photo 
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emission data (XPS). In figure [2l we have plotted the 
electronic density of states along with the experimen- 
tal data. The electronic density of states of the ECMR 
model (shown in the inset) are in good agreement with 
the XPS data, and establishes further credibility of our 
method. 

We have also applied our ECMR approach to model 
hydrogenated amorphous silicon. For this purpose, we 
started with a pure amorphous silicon configuration ob- 
tained via Reverse Monte Carlo simulation. The partic- 
ular RMC scheme to generate amoprhous silicon is dis- 
cussed in Ref. 13. The initial size of the model is 500 



atoms, which is then hydrogenated following a method 
similar to but not identical to Holender and Morgan . 
In our scheme, the dangling bonds are passivated by in- 
troducing H atoms in the RMC-generated continuous 
random network. Once the dangling bonds are identi- 
fied, passivation is done by placing a H atom at a dis- 
tance of 1.45 A to 1.65 A. The atoms are placed along 
the direction vector opposite to the sum of the three vec- 
tors connecting the central atom and its three neighbors. 
Once the passivation is done, the density of the model is 
adjusted to experimental density and the system is sub- 
jected to ECMR iteration to make the resulting config- 
uration consistent with both experimental data and the 
total energy. Any additional dangling bonds generated 
during ECMR iteration are also passivated, and the pro- 
cess is repeated until the dangling bonds cease to exist. 
The resulting model of hydrogenated amorphous silicon 
consists of 540 atoms with 7.4% H atom in the network. 
The model has defect concentration of about 2%, which 
consists of 5-fold coordinated floating bonds. 

The partial radial distribution functions for the model 
are plotted in figure [31 For the Si-Si case, most of the in- 
formation are already included in the starting model, and 
hence we do not expect much changes in the final model. 
For Si-H part, however, it is important to note that 
during hydrogenation and subsequent relaxation of the 
model via ECMR, the system does not change much of 
the pair correlations between Si and H atoms. The elec- 
tronic density of states (EDOS) for this 540-atom model 
is plotted in figure IH The EDOS is obtained using the 
first-principles density functional code Siesta within the 
local density approximation. The model clearly shows 
the presence of a gap in the spectrum along with a few 
tail states near the band edges. This is due to presence 
of few 5-fold floating bonds in the network. As a further 
test of our model, we have also computed the vibrational 
density of states (VDOS) within the harmonic approxi- 
mation by constructing the dynamical matrix using the 
electronic forces. The eigenvalues of the dynamical ma- 
trix are related to the square of the angular frequency 
of vibrations. The VDOS in figure H] distinctly show the 
presence of acoustic and optical peaks due to vibration 
of Si atoms. The high frequency vibrations in the den- 
sity of states correspond to the H atom movement in the 



network. 



CONCLUSION 

Wc present an inverse approach to construct atomistic 
models of materials using a combination of experimental 
data and a suitable first-principles force-field. Using a 
generalised penalty function, we merge both the power of 
density functional theory with experimental data so that 
the configurations are not only a minimum (metastable 
or global) of an energy functional but also consistent 
with available experimental data. The method is ap- 
plied to construct model configurations for two techno- 
logically important materials: glassy GeSe2 and hydro- 
genated amorphous silicon. The structural, electronic 
and vibrational properties of the models are compared 
with experiments, and have been observed to be in very 
good agreement with experimental data. 
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FIG. 1: Neutron weighted static structure factors for the 
ECMR, reverse Monte Carlo (RMC) and decorate-and-relax 
(PR) rnodel. The resuh for the DR model is taken from the 
Ref. Hol . The low wave vector part is magnified in the inset 
for comparison to our model with experimental data and the 
DR model. 
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FIG. 2: The experimental and computed electronic density 
of states (inset) of g-GeSe2. The results for models obtained 
from the ECMR, reverse Monte Carlo (RMC), and decorate- 
and-relax (DR) are shown in inset. The results for the DR 
model are taken from Ref. [13 for comparison. 
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FIG. 3: The partial pair correlation functions for Si-Si (upper 
panel) and and Si-H (lower panel) atoms for the 540-atom 
ECMR model described in the text. 




FIG. 4: The electronic density of states of the 540-atom 
model of a-Si:H obtained from ECMR method within the local 
density approximation. The band gap is clearly visible in 
the spectrum along with a few tail states originated from the 
floating bonds present in the network 
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FIG. 5: The vibrational density of the 540-atom model a- 
Si:H obtained within the harmonic approximation. The high 
frequency vibrations correspond to H atom movement in the 
network. The acoustic and optical peaks are clearly visible in 
the vibrational spectrum. 



